library(Seurat)
library(clusterProfiler)
library(tidyseurat)
library(tidySingleCellExperiment)
library(ggpubr)
library(ggrepel)
library(scMerge)
library(tidyverse)

# load seurat objects -----------
sobj <- read_rds('CRC-I/data/zy2020_tumor10x.rds')

sg_sobj <- read_rds('CRC-I/data/seekgene/crc-starsolo-annotated.rds')

sg_sobj |> count(orig.ident, Platform, genotype)

data("segList", package = "scMerge")

# pro-tumor genes in myeloid cells ------------
myeloid_sg <- sg_sobj |>
  filter(zhang2020_main == 'Myeloid cell')

mmp_list <- myeloid_sg |>
  rownames() |>
  str_subset('^MMP\\d+$')

tumor_prolifer <- c('IL6','TNF','EGF')

cathepsin_list <- myeloid_sg |> rownames() |> str_subset('^CTS[B|C|F|H|K|L|O|S|V|X|W]$')

angiogen_list <- c('VEGFA','CXCL8','SEMA4D','PROK2')

immusuppr_list <- c('TGFB1','TGFB2','IL10','CD40LG','CD274','XBP1','NOS2','TPH1')

sg_sobj |>
  DotPlot(features = immusuppr_list, group.by = 'genotype') +
  RotatedAxis()

myeloid_sg |>
  DotPlot(features = mmp_list, group.by = 'zhang2020_fine') +
  RotatedAxis()

myeloid_sg |>
  DotPlot(features = cathepsin_list, group.by = 'zhang2020_fine') +
  RotatedAxis()

myeloid_sg |>
  DotPlot(features = angiogen_list, group.by = 'zhang2020_fine') +
  RotatedAxis()

myeloid_sg |>
  DotPlot(features = immusuppr_list, group.by = 'zhang2020_fine') +
  RotatedAxis()

myeloid_sg |>
  VlnPlot(features = 'FCGR2B', group.by = 'zhang2020_fine', pt.size = 0)

sobj |>
  filter(Global_Cluster == 'Myeloid cell') |>
  VlnPlot(features = 'FCGR2B', group.by = 'Sub_Cluster', pt.size = 0)

## in latent clusters --------
myeloid_sg |>
  DotPlot(features = mmp_list, group.by = 'latent_cluster') +
  RotatedAxis()

myeloid_sg |>
  DotPlot(features = cathepsin_list, group.by = 'latent_cluster') +
  RotatedAxis()

myeloid_sg |>
  DotPlot(features = angiogen_list, group.by = 'latent_cluster') +
  RotatedAxis()

myeloid_sg |>
  DotPlot(features = immusuppr_list, group.by = 'latent_cluster') +
  RotatedAxis()
  

# scMerge all myeloid cells -----------
zh_tam <- sobj |>
  filter(str_detect(Global_Cluster, 'Myeloid cell'))

gene_merged <- tibble(gene = rownames(myeloid_sg)) |>
  filter(gene %in% rownames(zh_tam)) |>
  pull(gene)

# visualize CD45 expression as an ref to see if expression scale is the same between 2 data sets.
mrg_tam <- merge(zh_tam, myeloid_sg)

# keep only genes appeared in both datasets
mrg_tam <- subset(mrg_tam, features=gene_merged)

## semi-supervised scMerge ------
mrg_tam <- mrg_tam |> mutate(zhang2020_fine = ifelse(
  is.na(zhang2020_fine),
  str_remove(Sub_Cluster, '.+_'),
  zhang2020_fine))

system.time(scm_sce <- mrg_tam %>%
              as.SingleCellExperiment() %>%
              scMerge(ctl = segList$human$human_scSEG,
                      assay_name = 'scMerge',
                      hvg_exprs = "logcounts",
                      kmeansK = c(9, 9),
                      cell_type = mrg_tam$zhang2020_fine,
                      batch_name = 'Platform',
                      verbose = TRUE))

convert_neg2zero <- function(n){
  ifelse(n<0, 0, n)
}

sparsify_neg2zero <- function(mat){
  mat@x <- mat@x |> convert_neg2zero()
  
  mat
}

scm_sobj <- assay(scm_sce, 'scMerge') |>
  sparsify_neg2zero() |>
  CreateSeuratObject()

scm_sobj@meta.data <- mrg_tam@meta.data

up_seek <- scm_sobj |>
  FindMarkers(ident.1 = 'TT', group.by = 'genotype') |>
  as_tibble(rownames = 'gene')

up_seek |>
  filter(gene %in% mmp_list)

up_seek <- up_seek |>
  mutate(type = case_when(avg_log2FC > 1 & p_val_adj < .05 ~ 'Up in TT',
                          avg_log2FC < -1 & p_val_adj < .05 ~ 'Down in TT',
                          .default = 'NS'))

up_seek |>
  count(type)

up_seek_name <- up_seek |>
  filter(type != 'NS')

up_seek |>
  ggplot(aes(avg_log2FC, -log10(p_val_adj), color = type)) +
  geom_point() +
  geom_vline(xintercept = c(-1,1), linetype = 'dashed') +
  geom_hline(yintercept = 5, linetype = 'dashed') +
  geom_text_repel(data = up_seek_name, aes(label = gene), color = 'black') +
  scale_color_manual(values = c('blue','grey','red') ) +
  theme_pubr() +
  labs_pubr() +
  labs(title = 'Myeloid cells DEG: TT vs II+IT') +
  expand_limits(x = c(-6,6))

# test for SPP1+ TAM --------
Idents(scm_sobj) <- 'zhang2020_fine'

## TT vs other ------------
up_seek_spp1 <- scm_sobj |>
  FindMarkers(ident.1 = 'TT',
              group.by = 'genotype',
              subset.ident = 'TAM-SPP1') |>
  as_tibble(rownames = 'gene')

up_seek_spp1 |>
  filter(gene %in% mmp_list)

up_seek <- up_seek_spp1 |>
  mutate(type = case_when(avg_log2FC > 1 & p_val_adj < .05 ~ 'Up in TT',
                          avg_log2FC < -1 & p_val_adj < .05 ~ 'Down in TT',
                          .default = 'NS'),
         sets = case_when(gene %in% immusuppr_list ~ 'Immunosuppressive',
                          gene %in% c(cathepsin_list, mmp_list) ~ 'Tumor invasion',
                          gene %in% angiogen_list ~ 'Angiogenesis',
                          gene %in% tumor_prolifer ~ 'Tumor proliferation'
         ),
         alpha = case_when(is.na(sets) ~ .2, .default = 1))

up_seek |>
  count(type)

up_seek_name <- up_seek |>
  filter(type != 'NS')

up_seek_sets <- up_seek |>
  filter(!is.na(sets))

symmetry_x_lim <- ceiling(up_seek$avg_log2FC |> abs() |> max())

# preview 800*600 is ok
up_seek |>
  ggplot(aes(avg_log2FC, -log10(p_val_adj), color = type)) +
  geom_point() +
  geom_vline(xintercept = c(-1,1), linetype = 'dashed') +
  geom_hline(yintercept = 1.3, linetype = 'dashed') +
  geom_text_repel(data = up_seek_name, aes(label = gene), color = 'black') +
  scale_color_manual(values = c('blue','grey','red') ) +
  theme_pubr() +
  labs_pubr() +
  expand_limits(x = c(symmetry_x_lim, -symmetry_x_lim)) +
  labs(title = 'SPP1+ TAM DEG: TT vs II+IT')

up_seek |>
  ggplot(aes(avg_log2FC, -log10(p_val_adj), color = sets)) +
  geom_point(aes(alpha = alpha)) +
  geom_vline(xintercept = c(-.5,.5), linetype = 'dashed') +
  geom_hline(yintercept = 5, linetype = 'dashed') +
  geom_text_repel(data = up_seek_sets, aes(label = gene), color = 'black') +
  scale_color_manual(values = c('red','blue','orange','green2','grey') ) +
  expand_limits(x = c(symmetry_x_lim, -symmetry_x_lim)) +
  theme_pubr(legend = 'right', base_size = 22) +
  labs_pubr() +
  labs(title = 'SPP1+ TAM DEG: TT vs II+IT') +
  guides(alpha = 'none')

### GO enrich: upreg in TT -------------
go_res <- up_seek |>
  dplyr::filter(avg_log2FC > 0, p_val_adj < .05) |>
  pull(gene) |>
  enrichGO(OrgDb = 'org.Hs.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           minGSSize = 3,
           readable = TRUE)

go_res@result |>
  as_tibble() |>
  dplyr::filter(p.adjust < .05) |>
  slice_head(n = 10) |>
  mutate(Description = str_wrap(Description, width = 50) |>
           fct_reorder(Count)) |>
  ggplot(aes(Description, Count, fill = p.adjust)) +
  geom_col() +
  coord_flip() +
  theme_pubr(legend = 'right') +
  labs_pubr() +
  scale_fill_gradient(low = 'red', high = 'black') +
  labs(y = 'DEG count',
       title = 'GO BP pathways upregulated in TT vs II+IT SPP1+ TAM') +
  theme(plot.title.position = 'plot')

sig_geneid <- go_res@result |>
  as_tibble() |>
  filter(str_detect(Description, 'negative regulation of immune system process')) |>
  separate_longer_delim(geneID, '/') |>
  pull(geneID)


### downreg in TT -------------
go_res2 <- up_seek |>
  dplyr::filter(avg_log2FC < 0, p_val_adj < .05) |>
  pull(gene) |>
  enrichGO(OrgDb = 'org.Hs.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           minGSSize = 3,
           readable = TRUE)

go_res2@result |>
  as_tibble() |>
  dplyr::filter(p.adjust < .05) |>
  slice_head(n = 10) |>
  mutate(Description = str_wrap(Description, width = 50) |>
           fct_reorder(Count)) |>
  ggplot(aes(Description, Count, fill = p.adjust)) +
  geom_col() +
  coord_flip() +
  theme_pubr(legend = 'right') +
  labs_pubr() +
  scale_fill_gradient(low = 'blue', high = 'black') +
  labs(y = 'DEG count',
       title = 'GO BP pathways downregulated in TT vs II+IT SPP1+ TAM')+
  theme(plot.title.position = 'plot')

sig_geneid2 <- go_res2@result |>
  as_tibble() |>
  filter(Description == 'antigen processing and presentation') |>
  separate_longer_delim(geneID, '/') |>
  pull(geneID)

scm_spp <- scm_sobj |>
  dplyr::filter(str_detect(zhang2020_fine, 'SPP1')) |>
  ScaleData()

scm_spp |>
  DoHeatmap(features = c(sig_geneid,
                         sig_geneid2),
            group.by = 'genotype',
            group.colors = c('green3','blue','red'))

scm_spp |>
  DotPlot(features = list('up'=sig_geneid,'down'=sig_geneid2),
          group.by = 'Platform') +
  RotatedAxis()

up_seek |> dplyr::filter(gene %in% c('IL1B','IL6','TNF','TGFB1'))

scm_sobj |>
  DoHeatmap(features =  c(head(sig_geneid, n = 25),head(sig_geneid2, n=25)),
            group.by = 'genotype',
            group.colors = c('green3','blue','red'))

# scMerge all T cells -----------
zh_t <- sobj |>
  filter(str_detect(Global_Cluster, 'CD'))

gene_merged <- tibble(gene = rownames(myeloid_sg)) |>
  filter(gene %in% rownames(zh_tam)) |>
  pull(gene)

# visualize CD45 expression as an ref to see if expression scale is the same between 2 data sets.
mrg_tam <- merge(zh_tam, myeloid_sg)

# keep only genes appeared in both datasets
mrg_tam <- subset(mrg_tam, features=gene_merged)

## semi-supervised scMerge ------
mrg_tam <- mrg_tam |> mutate(zhang2020_fine = ifelse(
  is.na(zhang2020_fine),
  str_remove(Sub_Cluster, '.+_'),
  zhang2020_fine))

system.time(scm_sce <- mrg_tam %>%
              as.SingleCellExperiment() %>%
              scMerge(ctl = segList$human$human_scSEG,
                      assay_name = 'scMerge',
                      hvg_exprs = "logcounts",
                      kmeansK = c(9, 9),
                      cell_type = mrg_tam$zhang2020_fine,
                      batch_name = 'Platform',
                      verbose = TRUE))

scm_sobj <- assay(scm_sce, 'scMerge') |>
  CreateSeuratObject()

scm_sobj@meta.data <- mrg_tam@meta.data


# any other myeloid cell ------------
celltype_name <- 'CD1C+ cDC2'

celltype_id <- unique(scm_sobj$zhang2020_fine) |>
  str_subset('cDC2')

up_seek_spp1 <- scm_sobj |>
  FindMarkers(ident.1 = 'TT',
              group.by = 'genotype',
              subset.ident = celltype_id) |>
  as_tibble(rownames = 'gene')

## TT vs other ------------
up_seek <- up_seek_spp1 |>
  mutate(type = case_when(avg_log2FC > 1 & p_val_adj < .05 ~ 'Up in TT',
                          avg_log2FC < -1 & p_val_adj < .05 ~ 'Down in TT',
                          .default = 'NS'),
         sets = case_when(gene %in% immusuppr_list ~ 'Immunosuppressive',
                          gene %in% c(cathepsin_list, mmp_list) ~ 'Tumor invasion',
                          gene %in% angiogen_list ~ 'Angiogenesis',
                          gene %in% tumor_prolifer ~ 'Tumor proliferation'
         ),
         alpha = case_when(is.na(sets) ~ .2, .default = 1))

up_seek |>
  count(type)

up_seek_name <- up_seek |>
  filter(type != 'NS')

up_seek_sets <- up_seek |>
  filter(!is.na(sets))

symmetry_x_lim <- up_seek$avg_log2FC |> abs() |>
  max() |> ceiling()

## plot ------
up_seek |>
  ggplot(aes(avg_log2FC, -log10(p_val_adj), color = type)) +
  geom_point() +
  geom_vline(xintercept = c(-1,1), linetype = 'dashed') +
  geom_hline(yintercept = 2, linetype = 'dashed') +
  geom_text_repel(data = up_seek_name, aes(label = gene), color = 'black') +
  scale_color_manual(values = c('blue','grey','red'),
                     label = c('Down in TT','NS','Up in TT')) +
  expand_limits(x = c(-symmetry_x_lim,symmetry_x_lim)) +
  theme_pubr() +
  labs_pubr() +
  labs(title = str_glue(celltype_name,' DEG: TT vs II+IT'))

up_seek |>
  ggplot(aes(avg_log2FC, -log10(p_val_adj), color = sets)) +
  geom_point(aes(alpha = alpha)) +
  geom_vline(xintercept = c(-1,1), linetype = 'dashed') +
  geom_hline(yintercept = 5, linetype = 'dashed') +
  geom_text_repel(data = up_seek_sets, aes(label = gene), color = 'black') +
  scale_color_manual(values = c('red','blue','orange','green2','grey') ) +
  expand_limits() +
  theme_pubr(legend = 'right', base_size = 22) +
  labs_pubr() +
  labs(title = str_glue(celltype_name,' DEG: TT vs II+IT')) +
  guides(alpha = 'none')

### upreg in TT -------------
go_res <- up_seek |>
  dplyr::filter(avg_log2FC > 0.15, p_val_adj < .05) |>
  pull(gene) |>
  enrichGO(OrgDb = 'org.Hs.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           minGSSize = 3,
           qvalueCutoff = .05,
           readable = TRUE)

go_res@result |>
  as_tibble() |>
  dplyr::filter(p.adjust < .05) |>
  slice_head(n = 10) |>
  mutate(Description = str_wrap(Description, width = 42) |>
           fct_reorder(Count)) |>
  ggplot(aes(Description, Count, fill = p.adjust)) +
  geom_col() +
  coord_flip() +
  theme_pubr(legend = 'right') +
  labs_pubr() +
  scale_fill_gradient(low = 'red', high = 'black') +
  labs(y = 'DEG count',
       title = str_glue('GO BP pathways upregulated in TT vs II+IT ', celltype_name)) +
  theme(plot.title.position = 'plot')

sig_geneid <- go_res@result |>
  as_tibble() |>
  head(n = 1) |>
  separate_longer_delim(geneID, '/') |>
  pull(geneID)

### downreg in TT -------------
go_res2 <- up_seek |>
  dplyr::filter(avg_log2FC < -0.25, p_val_adj < .05) |>
  pull(gene) |>
  enrichGO(OrgDb = 'org.Hs.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           minGSSize = 3,
           qvalueCutoff = .05,
           readable = TRUE)

go_res2@result |>
  as_tibble() |>
  dplyr::filter(p.adjust < .05) |>
  slice_head(n = 10) |>
  mutate(Description = str_wrap(Description, width = 42) |>
           fct_reorder(Count)) |>
  ggplot(aes(Description, Count, fill = p.adjust)) +
  geom_col() +
  coord_flip() +
  theme_pubr(legend = 'right') +
  labs_pubr() +
  scale_fill_gradient(low = 'blue', high = 'black') +
  labs(y = 'DEG count',
       title = str_glue('GO BP pathways downregulated in TT vs II+IT ', celltype_name)) +
  theme(plot.title.position = 'plot')

sig_geneid2 <- go_res2@result |>
  as_tibble() |>
  filter(Description == 'antigen processing and presentation') |>
  separate_longer_delim(geneID, '/') |>
  pull(geneID)

scm_sobj |>
  DoHeatmap(features = c(sig_geneid,sig_geneid2),
            group.by = 'genotype',
            group.colors = c('green3','blue','red')) &
  theme(axis.text.y = element_blank())

up_seek |> dplyr::filter(gene %in% c('IL1B','IL6','TNF','TGFB1'))

scm_sobj |>
  DoHeatmap(features =  c('IL1B','IL6','TNF','TGFB1'),
            group.by = 'genotype',
            group.colors = c('green3','blue','red'),
            lines.width = 10)

